using System;
using System.Diagnostics;
using System.Collections.Generic;
using PortfolioAllocator.Optimization;

namespace PortfolioAllocator.Tests
{
    using NUnit.Framework;

    [TestFixture]
    class TestQPSolver
    {
        [Test]
        [Category("Calculation")]
        public void test1()
        {
            // int n = 2;   // size of x and matrix G
            // int p = 1;   // number of equalities
            // int m = 3;   // number of inequalities
            double[,] G = { { 4, -2 }, { -2, 4 } };
            double[] g0 = { 6, 0 };
            double[] x;

            double[,] CE = { { 1 }, { 1 } };
            double[] ce0 = { -3 };

            double[,] CI = { { 1, 1, 1 }, { 0, 0, 1 } };
            double[] ci0 = { 0, 0, 0 };

            QPSolver solver = new QPSolver();
            double min = solver.solve(G, g0, CE, ce0, CI, ci0, out x);

            double precision = 1e-10;
            Assert.AreEqual(12.00, min, precision, "Test1: min=12 FAILED");
            Assert.AreEqual(1, x[0], precision, "Test1: x0=1 FAILED");
            Assert.AreEqual(2, x[1], precision, "Test1: x1=2 FAILED");
        }

        [Test]
        [Category("Calculation")]
        public void test2()
        {
            int n = 3;   // size of x and matrix G
            int p = 0;   // number of equalities
            //int m = 3;   // number of inequalities
            double[,] G = { {1,0,0},
                            {0,1,0},
                            {0,0,1}};
            double[] g0 = { 0, -5, 0 };
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce0 = new double[p];

            double[,] CI = {{-4,2,0},
                            {-3,1,-2},
                            {0,0,1}};
            double[] ci0 = { 8, -2, 0 };

            QPSolver solver = new QPSolver();
            double min = solver.solve(G, g0, CE, ce0, CI, ci0, out x);

            double precision = 1e-10;
            Assert.AreEqual(-2.3809523809523809, min, precision, "test2: min=-2.38 FAILED");

            Assert.AreEqual(0.47619047619047616, x[0], precision, "test2: x0=0.47 FAILED");

            Assert.AreEqual(1.0476190476190477, x[1], precision, "test2: x1=1.04 FAILED");

            Assert.AreEqual(2.0952380952380953, x[2], precision, "test2: x2=2.09 FAILED");

        }

        [Test]
        [Category("Calculation")]
        public void test3()
        {
            //int n = 3;   // size of x and matrix G
            //int p = 1;   // number of equalities
            //int m = 3;   // number of inequalities
            double[,] G = {{56, 17, 31},
                   {17, 56, 11},
                   {31, 11, 23}};

            double[] g = { 0, 2, 1 };
            double[] x;

            double[,] CE = {{0},
                    {1},
                    {1}};
            double[] ce = { -1 };

            double[,] CI = {   {1,0,0},
                       {0,1,0},
                       {0,0,1}};
            double[] ci = { 0, 0, 0 };

            QPSolver solver = new QPSolver();
            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            /*! 
            \attention results from R: min= 9.017544 x = ( -0.10004247  0.03427774  0.16192441 )
            */
            double precision = 1e-10;
            Assert.AreEqual(11.438596491228068, min, precision, "min=11.43");
            Assert.AreEqual(0, x[0], precision, "x0=0");
            Assert.AreEqual(0.8070175438596493, x[2], precision, "x2=0.80");
        }

        [Test]
        [Category("Calculation")]
        public void test10x10Matrix()
        {
            int n = 10;   // size of x and matrix G
            int p = 1;   // number of equalities
            int m = 10;   // number of inequalities
            double[,] G = { // new double[n,n]
		                    {490, 525, 419, 391, 543, 560, 347, 433, 403, 447},
		                    {525, 689, 567, 554, 701, 764, 554, 748, 532, 755},
		                    {419, 567, 877, 548, 831, 781, 773, 830, 728, 1147},
		                    {391, 554, 548, 631, 704, 820, 602, 824, 659, 832},
		                    {543, 701, 831, 704, 949, 959, 817, 948, 774, 1089},
		                    {560, 764, 781, 820, 959,1107, 827,1098, 891, 1138},
		                    {347, 554, 773, 602, 817, 827, 848, 972, 707, 1116},
		                    {433, 748, 830, 824, 948,1098, 972,1355, 986, 1386},
		                    {403, 532, 728, 659, 774, 891, 707, 986,1061, 1172},
		                    {447, 755,1147, 832,1089,1138,1116,1386,1172, 1793}};
            double[] g = { 1, 2, 3, 4, 5, 1, 2, 3, 4, 5 };
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];
            // budget constraint
            CE[0, 0] = 1; CE[1, 0] = 1; CE[2, 0] = 1; CE[3, 0] = 1; CE[4, 0] = 1;
            CE[5, 0] = 1; CE[6, 0] = 1; CE[7, 0] = 1; CE[8, 0] = 1; CE[9, 0] = 1;
            ce[0] = -1;
            double[,] CI = new double[n, m];
            double[] ci = new double[m];

            double[,] tmp1 = {
		                    {1,0,0,0,0,0,0,0,0,0},
		                    {0,1,0,0,0,0,0,0,0,0},
		                    {0,0,1,0,0,0,0,0,0,0},
		                    {0,0,0,1,0,0,0,0,0,0},
		                    {0,0,0,0,1,0,0,0,0,0},
		                    {0,0,0,0,0,1,0,0,0,0},
		                    {0,0,0,0,0,0,1,0,0,0},
		                    {0,0,0,0,0,0,0,1,0,0},
		                    {0,0,0,0,0,0,0,0,1,0},
		                    {0,0,0,0,0,0,0,0,0,1},
                        	};
            for (int j = 0; j < m; ++j)
            {
                ci[j] = 0;
                for (int i = 0; i < n; ++i)
                    CI[i, j] = tmp1[i, j];
            }

            QPSolver solver = new QPSolver();
            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            double precision = 1e-10;
            Assert.AreEqual(229.22987096774222, min, precision, "min=228.22");

            double[] xR = {0.72354838709677616, -4.7622368429219653e-016, -1.0788462274402413e-017, 0.12425806451613025, 1.1686526676556317e-015, 
					 2.4798619458298718e-016, 0.15219354838709448, -5.4034334766373958e-018, 0, -1.0529191902618237e-017};
            for (int i = 0; i < n; ++i)
                Assert.AreEqual(xR[i], x[i], precision, "x[i=" + i + "] differs");

        }

        [Test]
        [Category("Calculation")]
        public void test20x20MatrixA()
        {
            int n = 20;   // size of x and matrix G
            int p = 1;   // number of equalities
            int m = 20;   // number of inequalities
            double[,] G = {
		            {  984, 1061, 825, 816,1126,1194, 954, 882, 886,  858,  301,  753, 1118,  679,  765,  592,  649,  861,  781, 1013},
		            { 1061, 1408,1142,1165,1479,1668,1386,1530,1288, 1492,  379,  913, 1278,  973, 1039,  994, 1007, 1289, 1147, 1462},
		            {  825, 1142,1706,1127,1695,1595,1609,1624,1858, 2216,  423,  753, 1012, 1108, 1301, 1385, 1461, 1706, 1599, 1707},
		            {  816, 1165,1127,1446,1568,1811,1367,1700,1527, 1634,  423,  957, 1157, 1066, 1133, 1261, 1217, 1450,  899, 1190},
		            { 1126, 1479,1695,1568,2090,2138,1882,1974,1971, 2155,  523, 1108, 1514, 1412, 1623, 1657, 1677, 1970, 1515, 1746},
		            { 1194, 1668,1595,1811,2138,2553,1916,2379,2130, 2285,  546, 1257, 1627, 1493, 1558, 1662, 1528, 1739, 1202, 1623},
		            {  954, 1386,1609,1367,1882,1916,2038,2041,1854, 2240,  503, 1013, 1365, 1291, 1521, 1545, 1474, 1804, 1662, 1757},
		            {  882, 1530,1624,1700,1974,2379,2041,2737,2314, 2745,  602, 1197, 1409, 1531, 1510, 1783, 1592, 1801, 1576, 1877},
		            {  886, 1288,1858,1527,1971,2130,1854,2314,2842, 2884,  629,  911, 1082, 1341, 1405, 1624, 1685, 1843, 1653, 1800},
		            {  858, 1492,2216,1634,2155,2285,2240,2745,2884, 3559,  674, 1065, 1193, 1495, 1582, 1916, 1928, 2203, 2222, 2332},
		            {  301,  379, 423, 423, 523, 546, 503, 602, 629,  674,  262,  422,  495,  422,  406,  437,  498,  540,  524,  671},
		            {  753,  913, 753, 957,1108,1257,1013,1197, 911, 1065,  422, 1051, 1236,  929,  886,  911,  840, 1019,  880, 1213},
		            { 1118, 1278,1012,1157,1514,1627,1365,1409,1082, 1193,  495, 1236, 1669, 1192, 1235, 1123, 1093, 1314, 1141, 1504},
		            {  679,  973,1108,1066,1412,1493,1291,1531,1341, 1495,  422,  929, 1192, 1305, 1282, 1405, 1376, 1487, 1208, 1525},
		            {  765, 1039,1301,1133,1623,1558,1521,1510,1405, 1582,  406,  886, 1235, 1282, 1462, 1503, 1498, 1692, 1336, 1487},
		            {  592,  994,1385,1261,1657,1662,1545,1783,1624, 1916,  437,  911, 1123, 1405, 1503, 1740, 1730, 1898, 1444, 1656},
		            {  649, 1007,1461,1217,1677,1528,1474,1592,1685, 1928,  498,  840, 1093, 1376, 1498, 1730, 2024, 2194, 1708, 1967},
		            {  861, 1289,1706,1450,1970,1739,1804,1801,1843, 2203,  540, 1019, 1314, 1487, 1692, 1898, 2194, 2593, 2076, 2237},
		            {  781, 1147,1599, 899,1515,1202,1662,1576,1653, 2222,  524,  880, 1141, 1208, 1336, 1444, 1708, 2076, 2356, 2361},
		            { 1013, 1462,1707,1190,1746,1623,1757,1877,1800, 2332,  671, 1213, 1504, 1525, 1487, 1656, 1967, 2237, 2361, 2935}
	            };
            double[] g = { 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0 };
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = { -1 };

            double[,] CI = new double[n, m];
            double[] ci = new double[m];
            for (int i = 0; i < n; ++i)
            {
                ci[i] = 0;
                CE[i, 0] = 1;   // budget constraint
                for (int j = 0; j < n; ++j)
                {
                    if (i == j)
                        CI[i, j] = 1;
                    else
                        CI[i, j] = 0;
                }
            }

            QPSolver solver = new QPSolver();

            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            double precision = 1e-10;
            Assert.AreEqual(133.05475500229471, min, precision, "min=133.05");
            double[] xR = {-6.9388939039072284e-018 , 4.1624865735667056e-018 , 0.000000e+00, -2.6145728905365803e-017, 1.9720655333901385e-015,
		            -2.456580041590864e-016, 1.1974118270962855e-017,7.635610461972103e-017, -7.5001233271516059e-017,9.1494695719120574e-017,
		            0.99973743892468592 , 1.4810700518439846e-016, -8.0036662075178434e-017 , 3.6312265982658726e-016, 0,
		            -2.9938807023611181e-017,  -4.620600402302014e-016,  2.3042054042459374e-016 , 6.2299765332939721e-017,  1.0720751191618808e-016};

            for (int i = 0; i < n; ++i)
                Assert.AreEqual(xR[i], x[i], precision, "x[i=" + i + "] differs");
        }

        [Test]
        [Category("Calculation")]
        public void test20x20MatrixB()
        {
            int n = 20;   // size of x and matrix G
            int p = 1;   // number of equalities
            int m = 21;   // number of inequalities
            double[,] G = {
		            {  984, 1061, 825, 816,1126,1194, 954, 882, 886,  858,  301,  753, 1118,  679,  765,  592,  649,  861,  781, 1013},
		            { 1061, 1408,1142,1165,1479,1668,1386,1530,1288, 1492,  379,  913, 1278,  973, 1039,  994, 1007, 1289, 1147, 1462},
		            {  825, 1142,1706,1127,1695,1595,1609,1624,1858, 2216,  423,  753, 1012, 1108, 1301, 1385, 1461, 1706, 1599, 1707},
		            {  816, 1165,1127,1446,1568,1811,1367,1700,1527, 1634,  423,  957, 1157, 1066, 1133, 1261, 1217, 1450,  899, 1190},
		            { 1126, 1479,1695,1568,2090,2138,1882,1974,1971, 2155,  523, 1108, 1514, 1412, 1623, 1657, 1677, 1970, 1515, 1746},
		            { 1194, 1668,1595,1811,2138,2553,1916,2379,2130, 2285,  546, 1257, 1627, 1493, 1558, 1662, 1528, 1739, 1202, 1623},
		            {  954, 1386,1609,1367,1882,1916,2038,2041,1854, 2240,  503, 1013, 1365, 1291, 1521, 1545, 1474, 1804, 1662, 1757},
		            {  882, 1530,1624,1700,1974,2379,2041,2737,2314, 2745,  602, 1197, 1409, 1531, 1510, 1783, 1592, 1801, 1576, 1877},
		            {  886, 1288,1858,1527,1971,2130,1854,2314,2842, 2884,  629,  911, 1082, 1341, 1405, 1624, 1685, 1843, 1653, 1800},
		            {  858, 1492,2216,1634,2155,2285,2240,2745,2884, 3559,  674, 1065, 1193, 1495, 1582, 1916, 1928, 2203, 2222, 2332},
		            {  301,  379, 423, 423, 523, 546, 503, 602, 629,  674,  262,  422,  495,  422,  406,  437,  498,  540,  524,  671},
		            {  753,  913, 753, 957,1108,1257,1013,1197, 911, 1065,  422, 1051, 1236,  929,  886,  911,  840, 1019,  880, 1213},
		            { 1118, 1278,1012,1157,1514,1627,1365,1409,1082, 1193,  495, 1236, 1669, 1192, 1235, 1123, 1093, 1314, 1141, 1504},
		            {  679,  973,1108,1066,1412,1493,1291,1531,1341, 1495,  422,  929, 1192, 1305, 1282, 1405, 1376, 1487, 1208, 1525},
		            {  765, 1039,1301,1133,1623,1558,1521,1510,1405, 1582,  406,  886, 1235, 1282, 1462, 1503, 1498, 1692, 1336, 1487},
		            {  592,  994,1385,1261,1657,1662,1545,1783,1624, 1916,  437,  911, 1123, 1405, 1503, 1740, 1730, 1898, 1444, 1656},
		            {  649, 1007,1461,1217,1677,1528,1474,1592,1685, 1928,  498,  840, 1093, 1376, 1498, 1730, 2024, 2194, 1708, 1967},
		            {  861, 1289,1706,1450,1970,1739,1804,1801,1843, 2203,  540, 1019, 1314, 1487, 1692, 1898, 2194, 2593, 2076, 2237},
		            {  781, 1147,1599, 899,1515,1202,1662,1576,1653, 2222,  524,  880, 1141, 1208, 1336, 1444, 1708, 2076, 2356, 2361},
		            { 1013, 1462,1707,1190,1746,1623,1757,1877,1800, 2332,  671, 1213, 1504, 1525, 1487, 1656, 1967, 2237, 2361, 2935}
	            };
            double[] g = { 1, 1, 1, 2, 2, 2, 3, 3, 3, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0 };
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];

            double[,] CI = new double[n, m];
            double[] ci = new double[m];

            for (int i = 0; i < n; ++i)
            {
                CE[i, 0] = 1; // budget constraint
                ci[i] = 0;
                for (int j = 0; j < n; ++j)
                {
                    if (i == j)
                        CI[i, j] = 1;
                    else
                        CI[i, j] = 0;

                }
                CI[i, n] = 0;
            }
            ce[0] = -1;
            CI[10, 20] = -1;  // x_10 < 0.5
            ci[20] = 0.5;

            QPSolver solver = new QPSolver();
            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            double precision = 1e-10;
            Assert.AreEqual(227.93962232006947, min, precision, "min=227");
            double[] xR = {0.39736676587258868 , 1.5955550630031551e-016 , -1.605970239081486e-017, -1.2714733230515607e-017, -2.3700270883563836e-015,
		            1.6326478502925091e-014, -8.1658493307301409e-017,-5.8875931146262893e-017, 1.8818504141217411e-017,-9.3417297014138661e-017,
		            0.50000000000000033 , 0.029773329717329104, -3.8328343392828524e-016 , 0 , -5.9875344402186207e-016,
		            0.072382776657123463, -3.2077850427122563e-016,  -6.0386673252785274e-017 , -3.9490204284560647e-018,  -9.0415599340273291e-018};

            for (int i = 0; i < n; ++i)
                Assert.AreEqual(xR[i], x[i], precision, "x[i=" + i + "] differs");
        }

        [Test]
        [Category("Calculation")]
        public void problemWithNoConstraints()
        {
            int n = 2;   // size of x and matrix G
            int p = 0;   // number of equalities
            int m = 0;   // number of inequalities
            double[,] G = { { 4, -2 }, { -2, 4 } };
            double[] g0 = { 6, 0 };
            double[] x;
            double[,] CE = new double[n, p];
            double[] ce0 = new double[p];

            double[,] CI = new double[n, m];
            double[] ci0 = new double[m];

            QPSolver solver = new QPSolver();

            double min = solver.solve(G, g0, CE, ce0, CI, ci0, out x);

            double precision = 1e-10;
            Assert.AreEqual(-6, min, precision, "min=-6");
            Assert.AreEqual(-2, x[0], precision, "x0=-2");
            Assert.AreEqual(-1, x[1], precision, "x1=-1");
        }

        [Test]
        [Category("Calculation")]
        public void testExceptionNonPositiveDefiniteMatrix()
        {
            int n = 2;   // size of x and matrix G
            int p = 2;   // number of equalities
            int m = 0;   // number of inequalities
            double[,] G = { {-1,1},
                            {1,1}};
            double[] g = { 0, 0 };
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];

            double[,] CI = new double[n, m];
            double[] ci = new double[m];

            QPSolver solver = new QPSolver();
            // Check: matrix is not positive definite exception
            Assert.Throws<System.Exception>(
                    delegate { double min = solver.solve(G, g, CE, ce, CI, ci, out x); },
                    "Trying ot solve problem with not-positive definite quadratic form");
        }

        [Test]
        [Category("Calculation")]
        public void testLinearlyDependentEqualityConstraints()
        {
            int n = 2;   // size of x and matrix G
            int p = 2;   // number of equalities
            int m = 0;   // number of inequalities
            double[,] G = { { 1, 0 }, { 0, 1 } };
            double[] g = new double[n];
            double[] x;

            double[,] CE = { { 1d, 1d }, { 1d, 1d } };   // Check: linearly dependent equality constraints
            double[] ce = new double[p];

            double[,] CI = new double[n, m];
            double[] ci = new double[m];

            QPSolver solver = new QPSolver();
            Assert.Throws<System.Exception>(
                delegate { double min = solver.solve(G, g, CE, ce, CI, ci, out x); },
                "Trying ot solve problem with Equality constraints are linearly dependent");
        }

        [Test]
        [Category("Calculation")]
        public void testInfeasibleProblem()
        {
            int n = 2;   // size of x and matrix G
            int p = 0;   // number of equalities
            //int m = 2;   // number of inequalities
            double[,] G = { { 1, 0 }, { 1, 1 } };
            double[] g = { 0, 1 };
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];

            double[,] CI = { { 1, 0 }, { 1, 0 } };
            double[] ci = { 1, -1 };

            QPSolver solver = new QPSolver();
            // Check: infeasible equality constraints
            Assert.Throws<System.Exception>(
                delegate { double min = solver.solve(G, g, CE, ce, CI, ci, out x); },
                "Trying ot solve problem with Infeasible problem");
        }

        [Test]
        [Category("Calculation")]
        public void noBindingConstraints()
        {
            int n = 1;   // size of x and matrix G
            int p = 0;   // number of equalities
            //int m = 2;   // number of inequalities

            //min (x-1)^2
            double[,] G = { { 2 } };
            double[] g = { -2 };
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];

            // subj to: x > 0 and x < 2
            double[,] CI = { { 1, -1 } };
            double[] ci = { 0, 2 };

            QPSolver solver = new QPSolver();
            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            double precision = 1e-10;
            Assert.AreEqual(-1, min, precision, "objective"); // should be 0.0 ????
            Assert.AreEqual(1, x[0], precision, "x0");
        }

        /// <summary>Small unconstrained test from: http://www.akiti.ca/QuadProgEx0Constr.html</summary>
        [Test]
        [Category("Calculation")]
        public void smallNoConstraints()
        {
            int n = 3;   // size of x and matrix G
            int p = 0;   // number of equalities
            int m = 0;
            double[,] G = { {5,-2,-1},
                    {-2,4,3},
                    {-1,3,5}};
            double[] g = { 2, -35, -47 };
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];

            double[,] CI = new double[n, m];
            double[] ci = new double[m];

            //min x^2
            QPSolver solver = new QPSolver();

            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            // There are no constraints
            double precision = 1e-10;
            Assert.AreEqual(-249, min, precision, "objective");
            Assert.AreEqual(3, x[0], precision, "x0");
            Assert.AreEqual(5, x[1], precision, "x1");
            Assert.AreEqual(7, x[2], precision, "x2");
        }

        /// <summary>Small test case from: http://www.me.utexas.edu/~jensen/ORMM/supplements/methods/nlpmethod/S2_quadratic.pdf</summary>
        [Test]
        [Category("Calculation")]
        public void smallWithConstraints()
        {
            int n = 2;   // size of x and matrix G
            int p = 0;   // number of equalities
            int m = 4;   // number of inequalities
            double[,] G = new double[n, n];
            double[] g = new double[n];
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];

            double[,] CI = new double[n, m];
            double[] ci = new double[m];

            G[0, 0] = 2; G[0, 1] = 0;
            G[1, 0] = 0; G[1, 1] = 8;
            g[0] = -8; g[1] = -16;

            CI[0, 0] = -1; CI[1, 0] = -1; ci[0] = 5;  // x1+x2 < 5
            CI[0, 1] = -1; CI[1, 1] = 0; ci[1] = 3;  // x1 < 3
            CI[0, 2] = 1; CI[1, 2] = 0; ci[2] = 0;  // x1 > 0
            CI[0, 3] = 0; CI[1, 3] = 1; ci[3] = 0;  // x2 > 0;

            QPSolver solver = new QPSolver();
            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            double precision = 1e-10;
            Assert.AreEqual(-31, min, precision, "objective");
            Assert.AreEqual(3, x[0], precision, "x0");
            Assert.AreEqual(2, x[1], precision, "x1");
        }

        [Test]
        [Category("Calculation")]
        public void dependentInequalityConstraints()
        {
            int n = 3;   // size of x and matrix G
            int p = 0;   // number of equalities
            int m = 8;   // number of inequalities
            double[,] G = new double[n, n];
            double[] g = new double[n];
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];

            double[,] CI = new double[n, m];
            double[] ci = new double[m];

            G[0, 0] = 5; G[0, 1] = -2; G[0, 2] = -1;
            G[1, 0] = -2; G[1, 1] = 4; G[1, 2] = 3;
            G[2, 0] = -1; G[2, 1] = 3; G[2, 2] = 5;
            g[0] = 2; g[1] = -35; g[2] = -47;

            CI[0, 0] = 1; CI[1, 0] = 0; CI[2, 0] = 0; ci[0] = -2; // x1 > 2
            CI[0, 1] = -1; CI[1, 1] = 0; CI[2, 1] = 0; ci[1] = 2; // x1 < 2
            CI[0, 2] = -1; CI[1, 2] = 0; CI[2, 2] = 0; ci[2] = 2; // x1 < 2
            CI[0, 3] = 1; CI[1, 3] = 0; CI[2, 3] = 0; ci[3] = -2; // x1 > 2

            CI[0, 4] = 0; CI[1, 4] = 1; CI[2, 4] = 0; ci[4] = -3; // x2 > 3
            CI[0, 5] = 0; CI[1, 5] = -1; CI[2, 5] = 0; ci[5] = 3; // x2 < 3
            CI[0, 6] = 0; CI[1, 6] = 1; CI[2, 6] = 0; ci[6] = -3; // x2 > 3
            CI[0, 7] = 0; CI[1, 7] = -1; CI[2, 7] = 0; ci[7] = 3; // x2 < 3

            QPSolver solver = new QPSolver();
            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            double precision = 1e-10;
            Assert.AreEqual(-245, min, precision, "objective");
            Assert.AreEqual(2, x[0], precision, "x0");
            Assert.AreEqual(3, x[1], precision, "x1");
            Assert.AreEqual(8, x[2], precision, "x2");
        }

        /// <summary>test on several constraints defining 2D extreme point </summary>
        [Test]
        [Category("Calculation")]
        public void degenerateExtremePoint()
        {
            int n = 2;   // size of x and matrix G
            int p = 0;   // number of equalities
            int m = 3;   // number of inequalities
            double[,] G = new double[n, n];
            double[] g = new double[n];
            double[] x;

            double[,] CE = new double[n, p];
            double[] ce = new double[p];

            double[,] CI = new double[n, m];
            double[] ci = new double[m];

            // f = x^2 + y^2
            G[0, 0] = 2; G[0, 1] = 0;
            G[1, 0] = 0; G[1, 1] = 2;
            g[0] = 0; g[1] = 0;

            CI[0, 0] = 1; CI[1, 0] = 0; ci[0] = -1;  // x > 1 
            CI[0, 1] = 0; CI[1, 1] = 1; ci[1] = -1;  // y > 1
            CI[0, 2] = 1; CI[1, 2] = 1; ci[2] = -2;  // x+y > 2

            QPSolver solver = new QPSolver();

            double min = solver.solve(G, g, CE, ce, CI, ci, out x);

            double precision = 1e-10;
            Assert.AreEqual(2, min, precision, "objective");
            Assert.AreEqual(1, x[0], precision, "x0");
            Assert.AreEqual(1, x[1], precision, "x1");
        }
    }

    [TestFixture]
    class TestAllocationSolver
    {
        [Test]
        [Category("Calculation")]
        public void findMinVariancePortfWeights()
        {
            double[,] cov = { {1.3 , 0.3}, {0.3, 1} };
            double[] r = { 0.15, 0.10 };
            double rMin = 0.13;
            double[] xMin = { 0, 0 };
            double[] xMax = { 1, 1 };
            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMinVariancePortfWeights(cov, r, rMin, xMin, xMax);

            double precision = 1e-10;
            Assert.AreEqual(1.0, x[0] + x[1], precision, "Position constraint");
            double expRet = aSolver.calcExpReturn(x, r);
            Assert.AreEqual(rMin, expRet, precision, "Expected return");
            double vol = aSolver.calcVolatility(x, cov);
            Assert.AreEqual(0.878635305459552, vol, precision, "Volatility");
            double lnRet = aSolver.calcExpLogReturn(x, r, cov);
            Assert.AreEqual(expRet - 0.5 * vol * vol, lnRet, precision, "Expected log return");

            double[] xExp = { 0.6, 0.4 };
            Assert.That(xExp, Is.EqualTo(x).Within(precision), "Min variance weights");
        }

        [Test]
        [Category("Calculation")]
        public void findMinVarianceLuenberger()
        {
            double[,] cov = {{2.3,	0.93,	0.62,	0.74,	-0.23},
                            {0.93,	1.4,	0.22,	0.56,	0.26},
                            {0.62,	0.22,	1.8,	0.78,	0.27},
                            {0.74,	0.56,	0.78,	3.4,	-0.56},
                            {-0.23,	0.26,	-0.27,	-0.56,	2.6}};
            double[] r = { 0.151, 0.125, 0.147, 0.0902, 0.1768 };
            double rMin = 0.14411780000;
            double[] xMin = { 0, 0, 0, 0, 0 };
            double[] xMax = { 1, 1, 1, 1, 1 };

            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMinVariancePortfWeights(cov, r, rMin, xMin, xMax);

            double precision = 1e-9;
            Assert.AreEqual(1.0, x[0] + x[1] + x[2] + x[3] + x[4], precision, "Position constraint");
            double expRet = aSolver.calcExpReturn(x, r);
            Assert.AreEqual(0.14411780000, expRet, precision, "Expected return");
            double vol = aSolver.calcVolatility(x, cov);
            Assert.AreEqual(0.816549341241325, vol, precision, "Volatility");
            double lnRet = aSolver.calcExpLogReturn(x, r, cov);
            Assert.AreEqual(expRet - 0.5 * vol * vol, lnRet, precision, "Expected log return");

            double[] xExp = { 0.136227796, 0.261905485, 0.227485167, 0.101867907, 0.272513645 };
            Assert.That(xExp, Is.EqualTo(x).Within(precision), "Min variance weights");
        }

        [Test]
        [Category("Calculation")]
        public void findMinVarianceSectorConstraints()
        {
            double[,] cov = {{2.3,	0.93,	0.62,	0.74,	-0.23},
                            {0.93,	1.4,	0.22,	0.56,	0.26},
                            {0.62,	0.22,	1.8,	0.78,	0.27},
                            {0.74,	0.56,	0.78,	3.4,	-0.56},
                            {-0.23,	0.26,	-0.27,	-0.56,	2.6}};
            List<PositionConstraint> constraints = new List<PositionConstraint>();
            PositionConstraint c;
            // 0 <= x[0] + x[1] <= 0.2
            c = new PositionConstraint(0, 0.20, new List<int>(new int[] {0,1}));
            constraints.Add(c);
            // 0.60 <= x[3] + x[4] <= 0.65
            c = new PositionConstraint(0.6, 0.65, new int[] { 3, 4 });
            constraints.Add(c);
            // 0.40 <= x[0] + x[2] + x[4] <= 0.5
            c = new PositionConstraint(0.4, 0.5, new List<int>(new int[] { 0, 2, 4 }));
            constraints.Add(c);
            // 0 <= x[0] <= 1
            c = new PositionConstraint(0, 1, new List<int>(new int[] { 0 }));
            constraints.Add(c);

            double precision = 1e-7;
            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMinVariancePortfWeights(cov, constraints);
            Assert.AreEqual(1.0, x[0] + x[1] + x[2] + x[3] + x[4], precision, "xA: Position constraint");

            double[] xExp = { 0, 0.2, 0.199740932642487, 0.3, 0.300259067357513 };

            Assert.GreaterOrEqual(x[0] + x[1], 0, "Constraint 1");
            Assert.LessOrEqual(x[0] + x[1], 0.20, "Constraint 2");
            Assert.GreaterOrEqual(x[3] + x[4], 0.60 - precision, "Constraint 3");
            Assert.LessOrEqual(x[3] + x[4], 0.65, "Constraint 4");
            Assert.GreaterOrEqual(x[0] + x[2] + x[4], 0.4 - precision, "Constraint 5");
            Assert.LessOrEqual(x[0] + x[2] + x[4], 0.5 + precision, "Constraint 6");
            Assert.GreaterOrEqual(x[0], 0, "Constraint 7");
            Assert.LessOrEqual(x[0], 1, "Constraint 8");
            Assert.That(x, Is.EqualTo(xExp).Within(precision), "Min variance weights");
        }

        [Test]
        [Category("Calculation")]
        public void findMinVarianceNoRetConstraint()
        {
            double[,] cov = {{2.3,	0.93,	0.62,	0.74,	-0.23},
                            {0.93,	1.4,	0.22,	0.56,	0.26},
                            {0.62,	0.22,	1.8,	0.78,	0.27},
                            {0.74,	0.56,	0.78,	3.4,	-0.56},
                            {-0.23,	0.26,	-0.27,	-0.56,	2.6}};
            double[] r = { 0.151, 0.125, 0.147, 0.0902, 0.1768 };
            double rMin = 0.0;
            double[] xMin = { 0, 0, 0, 0, 0 };
            double[] xMax = { 1, 1, 1, 1, 1 };

            double precision = 1e-9;
            AllocationSolver aSolver = new AllocationSolver();
            double[] xA = aSolver.findMinVariancePortfWeights(cov, r, rMin, xMin, xMax);
            Assert.AreEqual(1.0, xA[0] + xA[1] + xA[2] + xA[3] + xA[4], precision, "xA: Position constraint");
            double[] xB = aSolver.findMinVariancePortfWeights(cov, xMin, xMax);
            Assert.AreEqual(1.0, xB[0] + xB[1] + xB[2] + xB[3] + xB[4], precision, "xB: Position constraint");
            double[] xExp = { 0.109700102414048, 0.293619403874134, 0.222086143544066, 0.123027990564432, 0.25156635960332 };

            Assert.That(xA, Is.EqualTo(xB).Within(precision), "Min variance weights 1");
            Assert.That(xA, Is.EqualTo(xExp).Within(precision), "Min variance weights 2");
        }

        [Test]
        [Category("Calculation")]
        public void findMinvarianceUncorrelated()
        {
            double[,] cov = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
            double[] r = { 0.01, 0.02, 0.03 };
            double rMin = 0.0243d;
            double[] xMin = { 0, 0,  0 };
            double[] xMax = { 1, 1,  1 };

            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMinVariancePortfWeights(cov, r, rMin, xMin, xMax);

            double precision = 1e-9;
            Assert.AreEqual(1.0, x[0] + x[1] + x[2] , precision, "Position constraint");
            double expRet = aSolver.calcExpReturn(x, r);
            Assert.AreEqual(rMin, expRet, precision, "Expected return");
            double vol = aSolver.calcVolatility(x, cov);
            Assert.AreEqual(0.652520753182099, vol, precision, "Volatility");
            double lnRet = aSolver.calcExpLogReturn(x, r, cov);
            Assert.AreEqual(expRet - 0.5 * vol * vol, lnRet, precision, "Expected log return");

            double[] xExp = { 0.1183333333333, 0.3333333333333, 0.5483333333333 };
            Assert.That(x, Is.EqualTo(xExp).Within(precision), "Min variance weights");
        }

        [Test]
        [Category("Calculation")]
        public void findMinvarianceWithRiskFree()
        {
            double closeToZero = 1e-20;
            double[,] cov = {   {closeToZero,0,0},
                                {0,1,0},
                                {0,0,1}};
            double rRiskFree = 0.1;
            double[] r = { rRiskFree, 0.14, 0.2 };
            double rMin = 0;
            double[] xMin = { 0, 0, 0 };
            double[] xMax = { 1, 1, 1 };

            // Solution with rMin constraint
            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMinVariancePortfWeights(cov, r, rMin, xMin, xMax);

            double precision = 1e-9;
            Assert.AreEqual(1.0, x[0] + x[1] + x[2], precision, "Position constraint");
            double expRet = aSolver.calcExpReturn(x, r);
            Assert.AreEqual(rRiskFree, expRet, precision, "Expected return");
            double vol = aSolver.calcVolatility(x, cov);
            Assert.AreEqual(0.0, vol, precision, "Volatility");
            double lnRet = aSolver.calcExpLogReturn(x, r, cov);
            Assert.AreEqual(rRiskFree, lnRet, precision, "Expected log return");
            double[] xExp = { 1, 0, 0 };
            Assert.That(x, Is.EqualTo(xExp).Within(precision), "Min variance weights");

            // Solution without rMin constraint
            x[0] = 0;
            Assert.AreNotEqual(xExp, x, "weights should be reset now");
            x = aSolver.findMinVariancePortfWeights(cov, xMin, xMax);
            Assert.That(x, Is.EqualTo(xExp).Within(precision), "Min variance weights");
        }

        [Test]
        [Category("Calculation")]
        public void findMaxLogOptimalFromLuenberger()
        {
            double closeToZero = 1e-10;
            double[,] cov = {   {closeToZero,0,0,0},
                                {0,0.09, 0.02, 0.01},
                                {0,0.02, 0.07,-0.01},
                                {0,0.01,-0.01,0.03}};
            double rRiskFree = 0.1;
            double[] r = { rRiskFree, 0.24, 0.2, 0.15 };
            double[] xMin = { -100, 0, 0 , 0};
            double[] xMax = { 1, 3, 3 , 3};

            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMaxLogReturnPortfWeights(cov, r, xMin, xMax);

            double precision = 1e-7;
            Assert.AreEqual(1.0, x[0] + x[1] + x[2] + x[3], precision, "Position constraint");
            double expRet = aSolver.calcExpReturn(x, r);
            Assert.AreEqual(0.47420381357846908d, expRet, precision, "Expected return");
            double vol = aSolver.calcVolatility(x, cov);
            Assert.AreEqual(0.61172201201483656d, vol, precision, "Volatility");
            double lnRet = aSolver.calcExpLogReturn(x, r, cov);
            Assert.AreEqual(0.28710190358672916d, lnRet, precision, "Expected log return");
            double[] xExp = { (1 - 1.05 - 1.38 - 1.78), 1.05, 1.382165, 1.77707 };
            Assert.That(x, Is.EqualTo(xExp).Within(1e-3), "Min variance weights");
        }

        [Test]
        [Category("Calculation")]
        public void findMaxLogOptimalWithRiskFreeFromLuenberger()
        {
            double[,] cov = {   {0.09, 0.02, 0.01},
                                {0.02, 0.07,-0.01},
                                {0.01,-0.01,0.03}};
            double rRiskFree = 0.1;
            double[] r = { 0.24, 0.2, 0.15 };
            double[] xMin = { 0, 0, 0 };
            double[] xMax = {3, 3, 3 };

            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMaxLogReturnPortfWeights(cov, r, rRiskFree, xMin, xMax);
            double xSum = 0;
            for (int i =0; i < x.Length; ++i)
                xSum += x[i];
            double xRiskFree = 1 - xSum;
            double precision = 1e-7;
            double expRet = aSolver.calcExpReturn(x, r) + xRiskFree * rRiskFree;
            Assert.AreEqual(0.47420381357846908d, expRet, precision, "Expected return");
            double vol = aSolver.calcVolatility(x, cov);
            Assert.AreEqual(0.61172201201483656d, vol, precision, "Volatility");
            double lnRet = aSolver.calcExpLogReturn(x, r, cov) + xRiskFree * rRiskFree;
            Assert.AreEqual(0.28710190358672916d, lnRet, precision, "Expected log return");
            double[] xExp = { 1.05, 1.382165, 1.77707 };
            Assert.That(x, Is.EqualTo(xExp).Within(1e-3), "Min variance weights");
        }

        [Test]
        [Category("Calculation")]
        public void findMaxLogOptimalUncorrelated()
        {
            int n = 10;
            double[,] cov = new double[n, n];
            double[] r = new double[n];
            double[] xMin = new double[n];
            double[] xMax = new double[n];
            double volInstr = 0.2;
            double mu = 0.17;
            for (int i = 0; i < n; ++i)
            {
                cov[i, i] = volInstr * volInstr;
                r[i] = mu;
                xMin[i] = 0;
                xMax[i] = 1;
            }

            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMaxLogReturnPortfWeights(cov, r, xMin, xMax);

            double totX = 0;
            for (int i = 0; i < x.Length; ++i)
                totX += x[i];

            double precision = 1e-10;
            Assert.AreEqual(1.0, totX , precision, "Position constraint");

            double expRet = aSolver.calcExpReturn(x, r);
            Assert.AreEqual(mu, expRet, precision, "Expected return");

            double vol = aSolver.calcVolatility(x, cov);
            Assert.AreEqual(volInstr/Math.Sqrt(n), vol, precision, "Volatility");

            double lnRet = aSolver.calcExpLogReturn(x, r, cov);
            Assert.AreEqual(0.168, lnRet, precision, "Expected log return");

            double[] xExp = { 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1 };
            Assert.That(x, Is.EqualTo(xExp).Within(precision), "Min variance weights");
        }

        [Test]
        [Category("Calculation")]
        public void findMaxLogOptimal1SecAndCash()
        {
            double[,] cov = { { 1e-10, 0 }, { 0, 0.04 } };
            double[] r = { 0.1, 0.17 };
            double[] xMin = { -1, -1 };
            double[] xMax = { 2, 2 };

            AllocationSolver aSolver = new AllocationSolver();
            double[] x = aSolver.findMaxLogReturnPortfWeights(cov, r, xMin, xMax);

            double totX = 0;
            for (int i = 0; i < x.Length; ++i)
                totX += x[i];

            double precision = 1e-8;
            Assert.AreEqual(1.0, totX, precision, "Position constraint");

            double expRet = aSolver.calcExpReturn(x, r);
            Assert.AreEqual(-0.075 + 1.75*0.17, expRet, precision, "Expected return");

            double vol = aSolver.calcVolatility(x, cov);
            Assert.AreEqual(0.35, vol, precision, "Volatility");

            double lnRet = aSolver.calcExpLogReturn(x, r, cov);
            Assert.AreEqual(0.16125, lnRet, precision, "Expected log return");

            double[] xExp = { -0.75, 1.75 };
            Assert.That(x, Is.EqualTo(xExp).Within(precision), "Min variance weights");
        }
    }
}
